home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / NRPAS13 / DIFEQ.PAS < prev    next >
Pascal/Delphi Source File  |  1991-04-29  |  2KB  |  69 lines
  1. PROCEDURE difeq(k,k1,k2,jsf,is1,isf: integer; indexv: glindex;
  2.       ne: integer; VAR s: glsarray; nsi,nsj: integer;
  3.       y: glyarray; nyj,nyk: integer);
  4. (* Programs using routine DIFEQ must define the types
  5. TYPE
  6.    glindex = ARRAY [1..nyj] OF integer;
  7.    glsarray = ARRAY [1..nsi,1..nsj] OF real;
  8.    glyarray = ARRAY [1..nyj,1..nyk] OF real;
  9. and also declare the global quantities
  10. CONST
  11.    m=41;
  12. VAR
  13.    h,c2,anorm: real;
  14.    mm,n: integer;
  15.    x: ARRAY [1..m] OF real;
  16. in the main routine.   *)
  17. VAR
  18.    temp2,temp: real;
  19. BEGIN
  20.    IF (k = k1) THEN BEGIN
  21.       IF (((n+mm) MOD 2) = 1) THEN BEGIN
  22.          s[3,3+indexv[1]] := 1.0;
  23.          s[3,3+indexv[2]] := 0.0;
  24.          s[3,3+indexv[3]] := 0.0;
  25.          s[3,jsf] := y[1,1]
  26.       END ELSE BEGIN
  27.          s[3,3+indexv[1]] := 0.0;
  28.          s[3,3+indexv[2]] := 1.0;
  29.          s[3,3+indexv[3]] := 0.0;
  30.          s[3,jsf] := y[2,1]
  31.       END
  32.    END ELSE IF (k > k2) THEN BEGIN
  33.       s[1,3+indexv[1]] := -(y[3,m]-c2)/(2.0*(mm+1.0));
  34.       s[1,3+indexv[2]] := 1.0;
  35.       s[1,3+indexv[3]] := -y[1,m]/(2.0*(mm+1.0));
  36.       s[1,jsf] := y[2,m]-(y[3,m]-c2)*y[1,m]/(2.0*(mm+1.0));
  37.       s[2,3+indexv[1]] := 1.0;
  38.       s[2,3+indexv[2]] := 0.0;
  39.       s[2,3+indexv[3]] := 0.0;
  40.       s[2,jsf] := y[1,m]-anorm
  41.    END ELSE BEGIN
  42.       s[1,indexv[1]] := -1.0;
  43.       s[1,indexv[2]] := -0.5*h;
  44.       s[1,indexv[3]] := 0.0;
  45.       s[1,3+indexv[1]] := 1.0;
  46.       s[1,3+indexv[2]] := -0.5*h;
  47.       s[1,3+indexv[3]] := 0.0;
  48.       temp := h/(1.0-sqr(x[k]+x[k-1])*0.25);
  49.       temp2 := 0.5*(y[3,k]+y[3,k-1])-c2*0.25*sqr(x[k]+x[k-1]);
  50.       s[2,indexv[1]] := temp*temp2*0.5;
  51.       s[2,indexv[2]] := -1.0-0.5*temp*(mm+1.0)*(x[k]+x[k-1]);
  52.       s[2,indexv[3]] := 0.25*temp*(y[1,k]+y[1,k-1]);
  53.       s[2,3+indexv[1]] := s[2,indexv[1]];
  54.       s[2,3+indexv[2]] := 2.0+s[2,indexv[2]];
  55.       s[2,3+indexv[3]] := s[2,indexv[3]];
  56.       s[3,indexv[1]] := 0.0;
  57.       s[3,indexv[2]] := 0.0;
  58.       s[3,indexv[3]] := -1.0;
  59.       s[3,3+indexv[1]] := 0.0;
  60.       s[3,3+indexv[2]] := 0.0;
  61.       s[3,3+indexv[3]] := 1.0;
  62.       s[1,jsf] := y[1,k]-y[1,k-1]-0.5*h*(y[2,k]+y[2,k-1]);
  63.       s[2,jsf] := y[2,k]-y[2,k-1]-temp*((x[k]+x[k-1])
  64.          *0.5*(mm+1.0)*(y[2,k]+y[2,k-1])-temp2
  65.          *0.5*(y[1,k]+y[1,k-1]));
  66.       s[3,jsf] := y[3,k]-y[3,k-1]
  67.    END
  68. END;
  69.